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Abstract 

We briefly review simulation schemes for the investigation of rare transitions and 
we resume the recently introduced Transition Interface Sampling, a method in 
which the computation of rate constants is recast into the computation of fluxes 
through interfaces dividing the reactant and product state. 
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1 Rare events 



In many complex systems of physical importance transitions take place between 
stable states separated by a high (free) energy barrier. Examples are isomer- 
izations in clusters, chemical reactions, protein folding and crystal nucleation. 
Molecular simulations techniques such as molecular dynamics (MD) in principle 
enable the computation of the reaction rate constants, the search for transition 
states and the exploration of reaction mechanisms. But since the rate constant 
of the transition depends exponentially on the activation barrier height, the 
expectation time of a transition can easily become orders of magnitude longer 
than the molecular timescale which is usually measured in femtoseconds. Hence, 
when using straightforward MD the study of these rare events is far beyond cur- 
rent computer capabilities. 

In the past decade, a number of methods has been proposed that aim to 
tackle this timescale problem from different points of view. One class of methods 
focusses on escaping the initial state without making assumptions on the final 
state. This can be achieved by, for instance, artificially increasing the frequency 
of the rare event in a controlled way. The methods of Voter and collaborators 
follow this approach: hyperdynamics 0G] aims at lowering the energy difference 
between the top of the barrier and the initial basin, the parallel replica method 
[Sj exploits the power of parallel processing to extend the molecular simulation 
time, and temperature-accelerated dynamics speeds up the event by raising 
the temperature. The idea of driving energy into the system to escape the basin 
of the energy minimum in which the system is initially prepared is also at the 
basis of conformational flooding |H] , the Laio-Parrinello method |S] , and the 
enhanced sampling of a given reaction coordinate |9j. Another possible route 
is to coarse-grain the molecular dynamics on the fly and explore the resulting 
free-energy landscape ^UJ- Several methods are devoted to the exploration 
of the full potential energy surface through all its minima and saddle points. 
Examples are eigenvector following |11M12| . the activation- relaxation technique 
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of Barkcma and Mousseau, |13| . the dimer method of Hcnkclmann and Jonsson 
[H| . the kinetic Monte Carlo (MC) approach ^JEHEh anc ^ * ne discrete path 
sampling of Wales [TSl[TU| . 

When the initial and final state are known it is possible to generate paths 
connecting the two in the form of a discretized chain of states. This is the basis 
of a second class, the so-called two point boundary methods. One option is to 
find a minimal energy path on the potential energy surface, as in the Nudged 
Elastic Band method of Jonsson and collaborators [201 1211 1221 and in the 
string method of E et al. |24| . or to find a true dynamical path by minimizing 
a suitably chosen action |25j . Another possibility is to use modified stochastic 
equations of motion that guide the system from the initial to the final state |2fi| . 

To summarize, one can say that each of the above methods is well suited 
in one specific subclass of systems, but they also have their specific drawbacks. 
Some are inefficient for high dimensional systems or only give structural infor- 
mation and neglect the dynamics. Others are designed to find only one tran- 
sition or make heavily use on assumptions or prior knowledge of the system. 
In complex systems at finite temperature, concepts like the minimum energy 
(or action) path or the lowest saddle point are not very useful. The reaction 
is rather described by a ensemble of paths. Similarly, one cannot speak of a 
particular transition state but only of an ensemble of transition states. 

An important quantity describing the kinetics of rare events is the rate con- 
stant. The traditional way to calculate rate constants in complex condensed 
matter system is by the reactive flux method based on Transition State Theory 
[271 1281 12*31 1301 \'M\ . This method consists of two stages. First, the free energy 
is computed as a function of selected degrees of freedom describing the reac- 
tion from the initial to the final state (the reaction coordinate) , for instance by 
biased sampling techniques |32l \',VM IMj . This step is complemented with the 
calculation of a dynamical transmission coefficient, by starting short trajecto- 
ries from the maximum of the free energy barrier. However, in complex systems 
the correct reaction coordinate can be exceedingly difficult to find. If the reac- 
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tion coordinate does not capture the molecular mechanism, the biased sampling 
methods will suffer from substantial hysteresis when following the system over 
the barrier. Moreover, even if the free energy profile is obtained correctly for 
this particular (but wrong) reaction coordinate, the corresponding transmission 
coefficient will be very low, making an accurate evaluation problematic. 

To overcome these problems, Chandler and collaborators developed the Tran- 
sition Path Sampling (TPS) method. This technique gathers a collection of true 
dynamical trajectories connecting the states without any a priori assumption 
of the reaction coordinate. From the ensemble of pathways rate constants can 
be calculated and reaction mechanisms can be extracted [351 1361 \T7\ 138) . The 
method can be combined with parallel tempering |39j and stochastic dynamics 
can be used for the case of diffusive barriers 001 QD ■ Successful applications of 
TPS are, among others, ion pair dissociation in water [421 I43| . alanine dipeptide 
in vacuum and in aqueous solution |44| . neutral |45| and protonatcd 46, 42] wa- 
ter clusters, autoionization in water 0Hjj and the folding of a polypeptide |49| . 
For a detailed review on TPS see Refs. [501 15 II 152] . Similar techniques by Elbcr 
and Olender [531 1541 155) and Doniach et al jSE] sample discretized stochastic 
pathways based on the Onsager-Machlup action. Finally, we mention the topo- 
logical method of Tanase-Nicola and Kurchan [HZ] in which they suggest to use 
TPS in combination with saddle point searching vector walkers. 

In this paper we focus on transition path sampling, and in particular on the 
Transition Interface Sampling (TIS) method, a recent improvement over the 
TPS rate constant calculation [HE]- We briefly review the TIS method and give 
an example of its application. Further details on the derivations, algorithms 
and applications can be found in Ref. 58 . 

2 Transition Interface Sampling 

Consider a dynamical system prepared in the initial state A. The state is stable 
in the sense that trajectories will stay in that state for a time long compared to 
the molecular time scale (e.g. vibrations). Eventually, the trajectories cross the 
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barrier and reach the final state B. If the barrier is sufficiently high the system 
shows exponential relaxation and the rate constant ItAB is well defined. In what 
follows we assume that we can compute the evolution of trajectories in phase 
space, without further specifying the details of the dynamics. For instance, 
we could use deterministic integration of Hamilton's equations of motion or 
stochastic dynamics generated by a Langevin equation. 

The starting point of TIS is the partitioning of phase space by interfaces. 
We define an order parameter X(x) as a function of the phase space point x 
(consisting of positions and momenta of all particles in the system) and the 
interfaces z = . . . n as the hypersurfaces {ir|A(a;) = We assume that the 
interfaces do not intersect, that A<j_i < Ai, and we describe the boundaries 
of state A and B by Ao and A„ respectively. In the same spirit of TPS, the 
TIS method has the non-trivial advantage that the order parameter X(x) does 
not have to be a properly chosen reaction coordinate capturing the essence of 
the dynamical mechanism. Instead, it is sufficient that this function is able to 
characterize the basins of attraction of the stable states (HOj • The basis of the 
TIS method is the microscopic equation for the rate constant 

kAB = (<I>ti,o) / (Iia) ■ (1) 

Here (4>i,o) is the effective positive flux from interface - i.e. from state A - 
through interface i, where effective means that recrossings are not being counted. 
Equivalently, it is an average over all phase points that are first crossings through 
interface i but belong to trajectories that originated in A. The denominator 
(Iia) is a normalization factor taking into account all the phase points for which 
the corresponding trajectories come directly from A without having visited B 
(naturally, this includes the entire region A, and most part of the basin of 
attraction of A). In a simulation the ratio in Eq. Q is computed by starting 
a MD simulation in state A and counting the number of effective crosses with 
interface n, i.e. the number of times it reaches state B, per time unit. Since 
the transitions are rare, Eq. JQ) is of no practical use in this form, because 
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the simulation would have to be long enough to see at least one spontaneous 
event. We can improve efficiency by relating the effective flux (</>i,o) through an 
interface i to the effective flux (^.^o) through an interface i — 1 closer to A: 

(<fn,o) = (<pi-i,o) ■ PA(Xi\Xi-i)- (2) 

Here, PA(Xi\\j) is the conditional probability that a trajectory coming from A 
crosses interface i provided that it has passed interface j. Iteratively substituting 
Eq. ((2J in Eq. we can write 

k *B = ^ n^CA^lA,) = ^Pa(KW,), (3) 

which is the key equation of TIS. Similar equations can be written for the 
backward rate constant Uba- 

We compute the first term in Eq. by starting a simulation in A and 
counting the number of crossings with interface 1 per unit time. A statistically 
accurate value can be obtained by choosing the interface close enough to the 
initial state. The second term P/i(A n |Ai) is computed using the factorization 
in Eq. ©. First, we generate an ensemble of paths that start in A, cross 
interface 1 and eventually return to A or, alternatively, continue to cross the 
next interface 2. The probability Pa(A2|Ai) is the ratio of the number of paths 
that reach interface 2 to the total number of sampled paths. Subsequently, we 
generate an ensemble of paths starting in A and crossing interface 2, from which 
the next term -Pa (A3 IA2) can De obtained, and so on until we reach interface 
n and thus state B. Each step consists of a sampling of paths starting in A 
and crossing interface i = 1 . . . n — 1. This procedure is similar in spirit to the 
umbrella sampling technique in which different windows are employed to obtain 
free energy profiles • Note that the final rate constant Uab is independent of 
the choice of the interfaces as long as the first and last one are inside the basin 
of attraction of the stable states A and B respectively. 

The path ensembles are generated by constructing a random walk in path 
space employing a MC algorithm. Fundamental for both the TPS and TIS 
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methods is the MC move that generates a new path from an existing one, the 
shooting move . An existing path between the initial and final state is stored 
as collection of discrete time-slices. In a shooting move, a time slice is chosen at 
random, the momenta of all particles are slightly changed, and a new trajectory 
is created by integrating the system forward and backward in time. The new 
path is accepted with a probability chosen such that detailed balance is obeyed. 
In contrast to TPS, the path length does not have a prefixed value |36| . but is 
allowed to vary because the integration of the equations of motion is stopped 
when one of the two interfaces of the corresponding ensemble is reached 158] . 

3 Discussion 

We tested the TIS method on the isomerization of a diatomic molecule im- 
mersed in a repulsive fluid, a simple model system used before to test the TPS 
algorithms [3Sj. We report some details and the results in Fig. ^ more details 
can be found in Ref |58| . The forward rate constant for the isomerization corre- 
sponds to an average transition time k^g — (3.6 ± 0.1)s in real units for argon, 
which is indeed many orders of magnitude beyond the MD time-step ~ 4/s. In 
Ref. [SS] we showed that TIS is at least a factor 2 more efficient than TPS when 
computing the rate constant at equivalent conditions and same final relative 
error. The efficiency increases to 5 or more in systems with many recrossings. 

The TIS algorithm makes use of the MC moves developed for TPS, and in 
this sense TIS can be considered an extension of TPS. TIS retains the good 
features, such as the independence of prior knowledge of a reaction coordinate 
|5U| . and improves on the weaker points, for example, minimizing computational 
effort by allowing a variable path length. The spirit behind the TIS methodology 
for computing rate constants, however, is different. The concept of a positive 
effective flux gives a faster convergence because only positive terms contribute 
to the rate. The implementation of the computer algorithm becomes easier and 
one can apply multidimensional or discrete interface order parameters. These 
advantages make TIS more efficient in terms of computational effort. Finally, in 
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some cases a better characterization of true reactive paths can be achieved and 
non-true recrossings can be avoided through a proper choice of the interfaces, 
for instance, imposing kinetic energy constraints |58) . 

The TIS method has been successfully applied to two realistic cases, the 
folding of a polypeptide and hydration of ethylene [SHI- In this last case 
the method was combined with quantum ab-initio MD simulations. A recent 
variation of the TIS method for diffusive systems exploits very efficiently the loss 
of long time scale correlation by using a recursive reformulation of the crossing 
probability and the sampling of much shorter paths [60j . These results show that 
TIS is capable of studying rare event processes in complex systems efficiently and 
encourage even more challenging applications, such as isomerization in clusters 
and crystal nucleation, on which we plan to report in the future. 
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Figure 1: The dimer inter-particle potential (dashed line) and the TIS crossing 
probability Pa (A, Ai) (solid line) as function of the order parameter A = r, the 
dimer inter-particle separation. State A at the first minimum of the double-well 
corresponds to a compact state of the molecule and state B at second minimum 
to an extended state. The number of 2D solvent particles with unit diameter is 
25, the density is 0.7, the total energy was fixed at E=25 (all in Lennard- Jones 
reduced units [23). We computed the flux ($1,0) / (h A ) = 0.1196 ± 0.0005, and 
the crossing probability function Pa{X\Xi) by matching the partial functions 
Pa(K+i I Ai), i = 1 ... 5 in five windows of path ensemble simulations. The inter- 
faces in between are indicated by vertical lines. The error on the data points 
is within symbol size. The smooth line joining the points was created using a 
finer grid of interfaces . The horizontal plateau when approaching state B 
at A, i= 6 — 1-74 is an expression of the commitment of the trajectories to the 
final stable state. Paths that cross r ~ 1.7 always reach eventually the final 
interface without going back to A. The value of the plateau equals Pa(A„, Ai). 
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